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Abstract.  The  modeling  of  micro  flows  through  a  3D  micro  nozzle  was  investigated  by  means  of  the  DSMC  and  NS 
techniques  in  order  to  verify  the  influence  of  the  computational  scheme  parameters  on  the  predicted  flow  quantities.  The 
investigated  low  Reynold's  number  flow  has  a  complex  structure  containing  subsonic  and  supersonic  parts.  Since  the 
influence  of  the  boundary  layer  on  the  flow  characteristics  is  significant,  a  strong  effort  was  made  to  predict  the  boundary 
layer  structure  as  exactly  as  possible. 


1.  INTRODUCTION 

MEMS-based  propulsion  devices  are  under  development  in  order  to  deliver  precise  impulses  for  spacecraft 
maneuvers  Ref.  [1],  The  advantages  of  these  propulsion  devices  include  lightweight  materials  and  the  ability  to 
provide  versatile  thrust  levels  along  with  the  potential  of  clustering  of  such  thruster  devices.  A  wide  range  of 
phenomena  which  are  specific  to  the  MEMS  functionality  has  been  considered  in  Ref.  [2]  including  DSMC 
investigations  of  the  nozzle  flow  taking  into  account  the  removal  of  energy  from  the  flow  by  heat  transfer  to  the 
nozzle  wall.  Several  conclusions  were  obtained  from  this  work.  In  particular  it  was  found  that  a  strong  influence  of 
the  three  dimensionality  on  the  predicted  flow  properties  takes  place  due  to  the  specifics  of  the  nozzle  geometry. 

In  the  current  work,  a  strong  effort  is  undertaken  in  order  to  study  two  particular  3D  cases  of  a  similar  MEMS 
flow.  The  main  purpose  of  the  work  is  to  obtain  an  independence  of  the  calculated  flow  quantities  from  the 
numerical  scheme  parameters  for  both  cases.  The  cases  are  similar  except  for  the  stagnation  pressure  levels.  Both 
DSMC  and  NS  numerical  schemes  were  used  for  low  pressure  case  and  NS  scheme  was  used  for  high  pressure  case. 


2.  NOZZLE  GEOMETRY  AND  FLOW  CONDITIONS 

The  geometry  of  the  micro  nozzle  considered  in  this  work  is  the  same  as  developed  by  NASA  Glenn  researchers. 
A  configuration  of  the  nozzle  can  be  seen  in  Fig  1.  where  a  schematic  of  the  3D  computational  domain  for  a  series 
of  Navier  Stokes  (NS)  calculations  is  presented.  The  internal  dimensions  of  the  nozzle  are  presented  in  the  Table  1. 
The  high-temperature  flow  inside  the  micro  truster  is  generated  by  a  laser-ignited  solid  mono-propellant 
decomposition.  Molecular  nitrogen,  the  major  product  of  such  a  decomposition,  is  considered  in  the  current  studies. 
Two  different  stagnation  pressures  of  0.1  and  5  atm  are  studied.  The  stagnation  temperature  remains  equal  to  2,000 
K  for  both  cases  and  the  wall  temperature  is  fixed  at  300  K. 

Since  the  nozzle  flow  is  viscous  and  for  a  stagnation  pressure  of  0.1  atm,  portions  of  the  flow  will  be  in  the 
transitional  regime,  a  kinetic  technique  must  be  employed  to  accurately  predict  the  device  parameters  in  this  case. 
On  the  other  hand,  current  research  activities  are  aimed  at  increasing  the  device  performance  by  increasing  the 
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stagnation  pressure  in  the  device  to  levels  on  the  order  of  1-5  atm.  It  is  reasonable  to  use  the  NS  technique  for  these 
higher  pressure  levels. 


TABLE  1.  Internal  Nozzle  Dimension 


length  (sub  and  supersonic  parts  ): 

4930  / u  m 

inlet  height: 

600  m 

throat  height: 

300  ju  m 

exit  height: 

150  fj  in 

thickness  (third  dimension): 

600  jJ  m 

3.  NUMERICAL  METHODS 

The  flow  parameters  are  calculated  by  a  3D  DSMC  approach  for  a  stagnation  pressure  of  0.1  atm  and  by  a  3D  NS 
solver  for  a  stagnation  pressure  level  of  0. 1  and  5  atm.  It  must  be  noted  that  the  calculations  of  the  flow  parameters  at 
a  pressure  level  of  0.1  atm  is  challenging  for  both  methods.  The  flow  is  sufficiently  dense  for  numerical  accuracy  to 
be  of  concern  in  a  DSMC  calculation  (the  Kn  number  measured  across  the  wide  (third  dimension)  of  the  nozzle  inlet 
is  0.01)  and  it  becomes  even  more  dense  in  the  boundary  layer  close  to  the  cold  wall.  On  the  other  hand  the  flow 
regime  is  also  close  to  the  limits  that  are  solvable  by  a  NS  method  and  probably  slip  boundary  conditions  should  be 
used. 

The  DSMC  method  is  a  statistical  approach  of  the  Boltzmann  equation  describing  the  discrete  nature  of  the  gas. 
The  method  was  proposed  by  G.A.  Bird  Ref.  [3]  and  has  become  the  most  valuable  method  for  the  modeling  and 
simulation  of  rarefied  gases.  The  SMILE  Ref.  [4]  software,  which  implements  majorant  frequency  scheme,  of  Ref. 
[5],  is  used  in  these  studies. 

Implementation  of  the  DSMC  method  requires  time  and  space  discretization  by  sufficiently  small  quantities  to 
insure  that  the  predicted  flow  parameters  correspond  to  the  solution  of  the  Boltzmann  equation.  Furthermore  these 
numerical  parameters  are  adaptive  to  the  flow  gradients  and  a  reasonable  number  of  computational  particles  must  be 
involved  in  the  calculation  in  order  to  obtain  a  thrue  solution  of  the  Boltzmann  equation.  The  investigation  of  the 
influence  of  numerical  parameters  of  the  method  on  the  obtained  flow  quantities  is  the  goal  of  this  work. 

The  numerical  solution  of  the  Reynolds-averaged  Navier-Stokes  equations  for  viscous  flow  is  obtained  by  use  of 
the  GASP  Ref.  [6]  software  which  utilizes  a  finite  volume  descritization  of  the  computational  domain.  A  no-slip 
boundary  condition  is  used  to  model  the  gas-surface  interaction.  Both  subsonic  and  supersonic  parts  of  the  nozzle 
flow  are  modeled  in  the  calculations.  The  nozzle  inlet  Mach  number  was  chosen  to  be  0.1  in  order  to  satisfy  GASP 
requirements. 

The  3D  structural  mesh  is  created  and  optimized  using  GRIDGEN  Ref.  [7]  software  in  order  to  accommodate  the 
flow  gradients.  Several  steps  of  the  mesh  refinement  are  performed.  The  entire  computational  domain  is  divided  into 
a  number  of  zones  and  different  meshing  strategies  are  used  for  each  zone.  Special  attention  is  paid  for  developing 
the  mesh  in  the  boundary  layer  and  nozzle  lip  zones  where  the  flow  gradients  are  expected  to  be  large.  The  mesh 
structure  has  been  adopted  following  the  flow  gradients  predicted  by  each  run  of  the  numerical  procedure.  Again  the 
goal  of  this  investigation  is  to  obtain  independence  of  the  predicted  flow  parameters  from  further  mesh  refinement. 


4.  COMPUTATIONAL  STRATEGY  AND  RESULTS 

The  investigations  of  the  0.1  atm  case  by  means  of  the  DSMC  method  and  5  atm  case  by  means  of  the  NS 
method  are  now  completed  with  final  results  obtained.  The  0.1  atm  NS  calculations  are  currently  underway, 
therefore,  the  presented  results  are  preliminary. 
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4.1  Direct  Simulation  Monte  Carlo  Results 


Several  calculations  have  been  done  with  the  most  intensive  ones  taking  about  two  weeks  to  perform  on  a  cluster 
of  46  Intel  Xenon  2.4  GHz  processors.  The  computational  time  step  has  been  chosen  to  be  the  smallest  possible  for 
the  available  machine  (single  precision)  so  that  its  influence  on  the  results  can  be  probably  excluded. 

The  first  DSMC  calculation  has  been  performed  with  a  reasonable  yet  relatively  small  amount  of  -1,300,000 
particles  and  -670,000  cells.  It  should  be  noted  that  SMILE  software  adopts  the  computational  grid  to  the  flow 
gradients  while  keeping,  at  the  same  time,  the  average  number  of  particles  of  about  3  to  4  in  a  computational  cell. 
Therefore,  further  grid  adaptation  is  only  possible  by  adding  more  particles  to  the  computational  domain.  Four 
computational  runs  with  a  gradual  increase  of  the  number  of  particles  were  performed.  The  numbers  of  particles  at 
each  run  were  1,300,000;  12,200,000;  65,200,000  and  130,000,000.  The  final  increase  is  only  a  factor  of  two  only 
since  the  limits  of  the  computational  resources  on  the  cluster  were  reached.  Figure  2  presents  the  flow  density 
contours  obtained  in  the  cases  A  and  D.  Table  2  summarizes  the  computational  scheme  parameters  for  all  of  the 
DSMC  calculations  presented  on  the  figures.  Even  though  the  general  shape  of  the  number  density  contours  is  the 
same  for  the  both  runs,  the  boundary  layer  contours  in  the  subsonic  part  of  the  nozzle  are  quite  different.  The 
discrepancy  between  the  simulation  runs  is  even  more  evident  in  the  temperature  and  flow  velocity  profiles. 

Figure  5  presents  translational  temperature  across  the  nozzle  in  the  middle  of  its  supersonic  part  and  Fig.  6 
presents  x-velocity  profiles  across  the  nozzle  exit.  Figures  7  and  8  present  translational  temperature  and  the  x 
component  velocity  along  the  nozzle  centerline  respectively.  Figures  6  and  8  also  shows  NS  results  obtained  by 
Case  H  for  purposes  of  comparison.  The  figures  5-8  suggest  that  a  better  grid  strategy  should  be  employed  in  future 
modeling  efforts.  The  following  conclusion  could  be  made  by  analyzing  the  Figs.  5-8.  Even  though  the  essential 
DSMC  requirements  for  the  time  step  and  the  number  of  particles  in  a  computational  cell  are  satisfied  in  all  the 
calculations,  and  all  the  computational  runs  reach  steady  state,  quantitative  and  even  qualitative  changes  of  the  flow 
parameters  are  evident  when  advancing  from  a  low-resolution  mesh  (Case  A)  to  a  higher-resolution  mesh  (Case  D). 


4.2  Navier-Stokes  results 

The  initial  mesh  of  the  presented  NS  runs  has  spatial  resolution  similar  to  the  grid  used  in  the  first  DSMC 
calculation  (case  A).  About  750,000  mesh  points  are  deployed  in  the  computational  domain  and  special  attention  is 
paid  to  the  meshing  of  the  boundary  layer  and  the  nozzle  lip  regions.  The  refined  mesh  structure  has  about 
5,300,000  mesh  points  Table  3  summarizes  the  computational  scheme  parameters  for  all  of  the  NS  calculations. 

Figure  3  shows  the  Mach  contours  in  the  XY  symmetry  plane  for  the  3D  NS  calculation  for  the  initial  conditions 
at  0.1  atm.  This  result  was  obtained  using  coarse  mesh  structure  (Case  G)  It  can  be  seen  in  the  figure  that  the 
calculated  flow  is  incorrect,  even  qualitatively.  The  boundary  layer  occupies  practically  the  entire  supersonic  part  of 
the  nozzle.  Figure  4  shows  distribution  of  the  same  parameter  obtained  using  a  dense  mesh  structure  (Case  H).  The 
quality  of  the  parameter  distribution  became  more  reasonable.  X-velocity  profiles  across  the  nozzle  exit  and  along 
the  nozzle  centerline  for  this  case  are  presented  on  Figs.  6  and  8  respectively  in  comparison  with  the  same 
parameters  obtained  by  DSMC  method. 

Figure  9-12  present  Mach  number  contours,  temperature  contours,  x-velocity  profiles  across  the  nozzle  exit  and 
temperature  along  the  nozzle  centerline  respectively  for  the  inlet  pressure  of  5.0  atm.  The  increase  of  the  stagnation 
pressure  makes  the  conditions  more  suitable  for  the  NS  solver  and  better  results  can  be  expected.  It  is  seen  in  the 
figures  9-10  that  the  contours  of  parameters  are  qualitative  correct. 

Furthermore  figures  11-12  show  independence  of  the  results  obtained  for  the  inlet  pressure  of  5.0  atm  from  the 
computational  mesh  and  from  the  order  of  the  scheme  accuracy.  Solution  of  the  first  order  of  accuracy  obtained  on 
the  coarse  mesh  structure  changes  only  slightly  along  with  refinement  of  the  mesh  and  improving  the  scheme  order 
of  accuracy.  Unfortunately,  a  comparison  with  the  DSMC  solution  for  these  initial  conditions  is  not  possible  since 
the  DSMC  scheme  is  prohibitively  computationally  expensive  for  such  high  pressures. 
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4.  CONCLUSIONS 


The  numerical  studies  of  the  low  Reynolds  MEMS  flow  indicate  that  the  modeling  of  such  flows  remains 
challenging.  According  to  the  flow  characteristics,  a  NS  solver  should  be  suitable  to  the  modeling  of  such  flows, 
however,  our  initial  calculations  demonstrate  that  the  application  of  NS  to  the  MEMS  flows  for  transitional  Kn 
numbers  is  not  straightforward.  The  DSMC  modeling  is  shown  to  be  possible,  but  the  computational  cost  remains 
high.  It  is  not  clear  at  this  point  which  method  is  more  applicable  for  modeling  of  this  type  of  flows.  Further  work,  in 
which  a  NS  mesh  with  a  more  accurate  representation  of  the  flow  gradients  will  be  implemented,  may  reveal  better 
applicability  of  the  NS  solver.  Slip  jump  boundary  conditions  used  in  NS  will  allow  extend  the  applicability  of  the 
method.  Several  technical  improvements  are  also  planned  to  make  the  DSMC  method  faster  and  save  computational 
cost. 
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TABLE  2.  DSMC  scheme  parameters  TABLE  3.  NS  scheme  parameters 


Case 

Number  of  particles 

Number  of  cells 

Case 

Mesh  points 

Order  of  accuracy 

Pressure  [ 

A 

1,300,000 

670,000 

E 

750,000 

1 

5.0 

B 

12,000,000 

5,600,000 

F 

5,300,000 

2 

5.0 

C 

65,000,000 

27,000,000 

G 

750,000 

1 

0.1 

D 

130,000,000 

52,000,000 

H 

5,300,000 

1 

0.1 
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FIGURE  1.  Computational  Domain  for  a  series  of  NS  FIGURE  2.  Number  density  contours  (1  =  5.0  e+23) 

calculations  top:  case  A.  bottom:  case  D 

0.004 

Mach  Number  Contours 


FIGURE  3.  Mach  number  contours.  Case  G  FIGURE  4.  Mach  number  contours,  Case  H 


y  [M]  y  [M] 

FIGURE  5.  Translational  temperature  across  the  nozzle  in  the  FIGURE  6.  X-velocity  profiles  across  the  nozzle  exit.  The 
middle  of  the  supersonic  part.  The  nozzle  axis  is  located  at  x=0.  nozzle  axis  is  located  at  x=0. 
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FIGURE  7.  Translational  temperature  along  the  nozzle 
centerline. 


FIGURE  8.  X-velocity  along  the  nozzle 
centerline. 

T  0.004 


Temperature  Contours  [K] 


FIGURE  9.  Mach  number  contours.  Case  F 


FIGURE  10.  Temperature  contours.  Case  F 


FIGURE  11.  X-velocity  profiles  across  the  nozzle  exit.  The 
nozzle  axis  is  located  at  x=0. 


FIGURE  12.  Temperature  along  the  nozzle  centerline. 
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